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A non-slip boundary condition at a wall for the lattice Boltzmann method is 
presented. In the present method unknown distribution functions at the wall are 
assumed to be an equilibrium distribution function with a counter slip velocity 
which is determined so that fluid velocity at the wall is equal to the wall velocity. 
Poiseuille flow and Couette flow are calculated with the nine-velocity model to 
demonstrate the accuracy of the present boundary condition. 
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Recently, the lattice Boltzmann (LB) method 1-4 has been used for many kinds of simulations 
of viscous flows. In particular, the LB method has been successfully applied to problems of fluid 
flows through porous media 5,6 and multiphase fluid flows. 7,8 On the other hand, as for the 
implementation of a non-slip boundary condition several approaches have been proposed. 9-14 
The bounce-back boundary condition 9 has been usually used to model stationary walls. However, 
it has been found that the bounce-back boundary condition has errors in velocity at the wall. 11-14 
Skordos 10 proposed a method for calculating particle distributions at a boundary node from fluid 
variables with the gradients of the fluid velocity. In his method the density is assumed to be 
known at the boundary. Noble et al. 11 developed a method for calculating the density at the 
boundary and the unknown components of the particle distributions. While their method gives 
accurate results with the seven-velocity model, it is not clear whether the method can be applied 
to other velocity models. 

In this Letter, a new approach for applying the non-slip boundary condition at the wall 
is presented. For explanation we use the LB method with the BGK collision model. 3,4 In the 
method the evolution of the distribution function /j(x, t) of particles with velocity at the 
point x and at time t is computed by the following equation: 

f i (x + *At,t + At) = -\[h(^t) - /Hx,*)], (1) 

where /^(x, t) is an equilibrium distribution function, r is a single relaxation time, and At is 
a time step during which the particles travel a grid spacing. The fluid mass density p and the 
fluid velocity u are defined in terms of the particle distribution function by 

p = E/*. ( 2 ) 

i 

"=-E/.c, (3) 

P i 

Here, we use the nine- velocity model 4,15 to explain the procedure of the method, but it is straight- 
forward to apply the method to other velocity models. The nine- velocity model has velocity vec- 
tors, c = 0, c{ = [cos(vr(i-l)/2),sin(7r(i-l)/2)], and cj 1 = y/2[cos(ir(i- \)/2), sin(7r(z- §)/2)] 
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for i = 1, ■ ■ ■ ,4. An equilibrium distribution function of this model is given by 
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At the non-slip wall we must specify the distribution functions of the particles whose velocity 
points to fluid region. In Fig. 1 the distribution functions f£, ff 1 , and f^ 1 are unknowns. In 
the kinetic theory of gases the assumption of diffuse reflection is often used at the wall. In 
this approximation gas molecules that strike the wall are assumed to leave it with a Maxwellian 
velocity distribution having the velocity and the temperature of the wall. In general, the velocity 
along the wall obtained with this assumption is not equal to that of the wall velocity, although 
the normal velocity to the wall is equal to that of the wall velocity. The difference between 
the fictitious velocity and the wall velocity is called the slip velocity. The idea of the present 
method is that the unknown distribution functions are assumed to be an equilibrium distribution 
function with a counter slip velocity which is determined so that the fluid velocity at the wall is 
equal to the wall velocity. That is, in the case of Fig. 1 the unknown distribution functions f£, 
fi 1 , and are assumed to be 
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where u w and i> w are the x and y components of the wall velocity, and p' and u' are unknown 
parameters. The unknown u' is the above-mentioned counter slip velocity. It is noted that we 
have no normal velocity jump at the wall because as mentioned above there exists no difference 
of the normal velocity to the wall between the fluid and the wall on the assumption of diffuse 
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reflection. The two unknown parameters are determined on the condition that the fluid velocity 
at the wall is equal to the wall velocity. Thus, we obtain two equations corresponding to the x 
and y components of the fluid velocity. Moreover, the density at the wall, p w , is an unknown 
quantity and is calculated by Eq. (2). Therefore, we finally obtain three equations for the three 
unknowns. The solutions are as follows: 

' 7o + // + /3 / + 2(/ 4 / + / 3 // + /i / )l, (io) 
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It is noted that the same procedure can be applied to other velocity models in spite of the 
number of particle velocities. In general, for isothermal models in D-dimensional space we have 
D+l constraints (D equations for each component of the velocity and one density equation at 
the wall) for D unknown parameters (p' and D— 1 slip velocities) and an unknown density at 
the wall. 

To demonstrate the accuracy of the present boundary condition, Poiseuille flow is calculated 
with the nine-velocity model. A two-dimensional steady flow between stationary parallel walls 
at y = —L and y = +L with a constant pressure gradient is considered. The pressure gradient 
is maintained by a density difference between inlet and outlet. At the inlet and the outlet the 
unknown distribution functions are determined as follows. 16 At the inlet, the unknown distribu- 
tion functions f{, ff 1 , and j\ l are calculated by adding a constant value to the corresponding 
distribution functions at the outlet so that f[\ m = //lout + C, // 7 |i n = /i 7 lout + \C, and 
// 7 1 in = /i 7 |out + \C with reference to the equilibrium distribution functions given by Eqs. (5) 
and (6). The constant value C is determined by setting the density calculated by Eq. (2) to be 
a given value at the inlet. The same procedure is used for calculating the unknown distribution 
functions at the outlet. In addition, the unknown distribution functions at four corners of the 
inlet and the outlet are calculated as follows. For example, at the lower corner of inlet, f{ and 
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fl 1 are calculated by the above-mentioned method. Then, f£, ff 1 , and fJ, 1 are calculated by 
the present non-slip boundary condition. The same method is used at the other three corners. 
In the following calculations 21 nodes are used between the walls. Figure 2 shows calculated 
velocity profiles u/« max (u mSuX is the velocity at y = 0) with the present boundary condition 
for various values of r. It is found that the results for 0.7 < r < 20 agree with the analytical 
solution within machine accuracy. The values of the counter slip velocity v! are shown in Fig. 
3. It is noted that the values of u' are negative for all r. The magnitude of the counter slip 
velocity increases as r becomes larger. Needless to mention, the magnitude of the counter slip 
velocity depends on the number of nodes between the walls. For comparison with the present 
results, the results with the bounce-back boundary condition in which j\ = /i,// 7 = fl 1 , and 
f] 1 = fl 1 are shown in Fig. 4. Also, the results under the condition with v! = in the present 
boundary condition, which corresponds to the usual diffuse reflection in the kinetic theory of 
gases, are shown in Fig. 5. It is seen that both results in Fig. 4 and Fig. 5 have slip velocities 
at the wall and the slip velocity increases as r becomes larger. In addition, we can see that the 
slip velocity under the condition with v! = in the present boundary condition is larger than 
that with the bounce-back boundary condition for the same value of r. This is because the 
bounce-back boundary condition has stronger non-slip effect than the usual diffuse reflection. 
In the next problem, we calculate Couette flow to demonstrate the accuracy of the LB method 
with the present boundary condition. In the problem, the upper plate at y = +L is moved with 
velocity U at t > 0, and the lower plate at y = — L is at rest. Figure 6 shows the calculated 
velocity profiles u/U at 200 time steps with 21 nodes between the walls and with r = 1. It is seen 
that the results agree well with the analytical solution. To determine the convergence rate, we 
perform simulations with 11, 21, 41, and 81 nodes between the walls. The errors are compared 
at the same dimensionless time corresponding to 200 time steps with 21 nodes. Figure 7 shows 
the error norms E\ = J2 y \ u ~ u *\/ J2 y \ u *\ an d E 2 = ^Jj2 y ( u ~ u*) 2 / \Jj2 y ( u *) 2 where u* is the 
analytical solution and the sums are taken over the same internal 9 nodes between the walls 
for all cases. The slope, m, of the convergence is m = 2.0004 for E±, and m = 2.0006 for E 2 . 
It is clearly found that the LB method with the present boundary condition is a second-order 
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scheme. 

From these results, we can conclude that the present boundary condition is accurate to model 
a non-slip flat boundary in the lattice Boltzmann simulations. It should be noted, however, that 
the present method has difficulties in dealing with corners. One of the possible ways to handle 
corners is to calculate the unknown distribution functions on each side at the corner and to 
average the results of common unknown distribution functions, while we should keep in mind that 
this approach produces errors due to rounded corners. More rigorous considerations for corners 
remain in future work. Although we use the nine- velocity model to explain the present boundary 
condition, the method can be directly applied to other velocity models; even to thermal models 
and to three-dimensional models. For thermal models, we may use an equilibrium distribution 
function with a counter temperature jump as well as the counter slip velocity. The value of the 
counter temperature jump can be determined by the condition for temperature at the wall. For 
three-dimensional models, we may use a counter slip velocity which has two components on the 
surface of the wall. The calculated results with the other velocity models will be reported in the 
future. 
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FIGURE CAPTION 



Fig. 1. Particle distribution functions of the nine- velocity model at the wall. 

Fig. 2. Calculated velocity profiles for Poiseuille flow with the present boundary condition for 
various values of r. 

Fig. 3. Counter slip velocity v! versus single relaxation time r in the calculations with the 
present boundary condition. 

Fig. 4. Calculated velocity profiles for Poiseuille flow with the bounce-back boundary condition 
for various values of r. 

Fig. 5. Calculated velocity profiles for Poiseuille flow under the condition with u' = in the 
present boundary condition for various values of r. 

Fig. 6. Calculated velocity profiles for Couette flow at 200 time steps with 21 nodes between 
the walls and with r = 1. 

Fig. 7. Error norms for Couette flow calculations with 11, 21, 41, and 81 nodes between the 
walls. 
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Fig. 1. Particle distrdfcaticn functions of the nine- 
velocity rrodel at the veil. 
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Fig. 2. Calculated velocity profiles for Fbiseuille 
flew with the present boundary condition 
for various values of 
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Fig. 3. Guitar slip velocity u' versus single reLascicn 
time in the calculations with the present bound- 
ary bounary condition. 
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Fig. 4. Calculated velocity profiles for Fbiseuille flow 
with the bounce-back boundary cordition for 
various values of 
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Fig. 5. Calculated velocity profiles for Fbiseullle flow 
under the oorriiticn with u' = in tine present 
bouncary oonditicn for various values of 
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Fig. 6. Calculated velocity profiles for Gouette flow 
at 200 tirre steps with 21 nodes between the 
walls and with = 1. 
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Fig. 7. Error nouns for Gcuette flow calculations with 
11, 21, 41, and 81 nodes between the walls. 
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